library(data.table)
library(magrittr)
library(parallel)
library(dplyr)
library(boot)
library(INLA)

setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM")
dir1="./bayes_pvalue_beta0/"
dir2="./bayes_BF/"
id=read.table("./20201120/at.least.one.AShM.in.DC.unitID.txt",head=T,sep="\t")
id=as.character(id$x)
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")
i=1
fn1=paste0(dir1,group1[i],".bayes_p.txt")
file1=read.table(fn1,head=T,sep = "\t")
file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
file1=data.frame(unitID=file1$unitID,normal_reads1=file1$normal_reads1,normal_reads2=file1$normal_reads2,tumor_reads1=file1$tumor_reads1,tumor_reads2=file1$tumor_reads2)
file1=file1[file1$unitID %in% id,]
names(file1)=c("unitID",paste(group1[i],names(file1)[2:5],sep="_"))
file=file1
for(i in 2:6){
  fn1=paste0(dir1,group1[i],".bayes_p.txt")
  file1=read.table(fn1,head=T,sep = "\t")
  file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
  file1=data.frame(unitID=file1$unitID,normal_reads1=file1$normal_reads1,normal_reads2=file1$normal_reads2,tumor_reads1=file1$tumor_reads1,tumor_reads2=file1$tumor_reads2)
  file1=file1[file1$unitID %in% id,]
  names(file1)=c("unitID",paste(group1[i],names(file1)[2:5],sep="_"))
  file=merge(file,file1,by="unitID",all=T)
}
filekp=file
refnum=grep(names(filekp),pattern="normal_reads1")###左边是ref，右边是alt
altnum=grep(names(filekp),pattern="normal_reads2")
file=filekp[,c(refnum,altnum)]	###control
result=matrix(,nrow(file),ncol=4)
for (i in 1:nrow(file)){
reads=as.numeric(as.matrix(file)[i,])
ref=reads[1:(ncol(file)/2)]
var=reads[(ncol(file)/2+1):ncol(file)]
tid=(1:6)-1
sl=which(!is.na(ref))
if ((length(sl)/2)>=1){
df=data.frame(y=var[sl],Totalreads=ref[sl]+var[sl],x2=tid[sl])
null <- logit(median(1-df[,1]/df[,2]))
formula = y ~ 1 + f(x2, model = "iid")
#formula = ALT_COUNT ~ 1 + f(TISSUE_ID, model = "iid")
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Totalreads,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}else{
df=data.frame(y=var[sl],Totalreads=ref[sl]+var[sl])
null <- logit(median(1-df[,1]/df[,2]))
formula = y ~ 1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Totalreads,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	#ci95 <- c(coef[4], coef[5])
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}
}
filekp$DC.con.pvalue=result[,1]
filekp$DC.con.FDR=p.adjust(filekp$DC.con.pvalue,method = "BH")
filekp$DC.con.beta0=result[,2]
filekp$DC.con.beta0.lower=result[,3]
filekp$DC.con.beta0.upper=result[,4]
rt=data.frame(unitID=filekp$unitID,filekp[,26:30])
write.table(rt,"./20201120/DC.con.beta0.txt",quote=F,row.names = F,sep = "\t")


filekp=file
refnum=grep(names(filekp),pattern="tumor_reads1")###左边是ref，右边是alt
altnum=grep(names(filekp),pattern="tumor_reads2")
file=filekp[,c(refnum,altnum)]	###case
result=matrix(,nrow(file),ncol=4)
for (i in 1:nrow(file)){
reads=as.numeric(as.matrix(file)[i,])
ref=reads[1:(ncol(file)/2)]
var=reads[(ncol(file)/2+1):ncol(file)]
tid=(1:6)-1
sl=which(!is.na(ref))
if ((length(sl)/2)>=1){
df=data.frame(y=var[sl],Totalreads=ref[sl]+var[sl],x2=tid[sl])
null <- logit(median(1-df[,1]/df[,2]))
formula = y ~ 1 + f(x2, model = "iid")
#formula = ALT_COUNT ~ 1 + f(TISSUE_ID, model = "iid")
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Totalreads,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}else{
df=data.frame(y=var[sl],Totalreads=ref[sl]+var[sl])
null <- logit(median(1-df[,1]/df[,2]))
formula = y ~ 1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Totalreads,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	#ci95 <- c(coef[4], coef[5])
result[i,1]=post_pred_p
result[i,2]=coef$mean
result[i,3]=coef$`0.025quant`
result[i,4]=coef$`0.975quant`
}
}
filekp$DC.case.pvalue=result[,1]
filekp$DC.case.FDR=p.adjust(filekp$DC.case.pvalue,method = "BH")
filekp$DC.case.beta0=result[,2]
filekp$DC.case.beta0.lower=result[,3]
filekp$DC.case.beta0.upper=result[,4]
rt=data.frame(unitID=filekp$unitID,filekp[,26:30])
write.table(rt,"./20201120/DC.case.beta0.txt",quote=F,row.names = F,sep = "\t")